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Abstract 



^ I This paper is concerned with SIR (susceptible -^ infected — )• removed) household 

epidemic models in which the infection response may be either mild or severe, with the 
type of response also affecting the infectiousness of an individual. Two different models 
are analysed. In the first model, the infection status of an individual is predetermined, 
perhaps due to partial immunity, and in the second, the infection status of an individual 
depends on the infection status of its infector and on whether the individual was infected 
by a within- or between-household contact. The first scenario may be modelled using 
^ ' a multitype household epidemic model, and the second scenario by a model we denote 

r^ ■ by the infector-dependent-severity household epidemic model. Large population results 

of the two models are derived, with the focus being on the distribution of the total 
^ numbers of mild and severe cases in a typical household, of any given size, in the event 

(^ ' that the epidemic becomes established. The aim of the paper is to investigate whether 

^D . it is possible to determine which of the two underlying explanations is causing the 

varying response when given final size household outbreak data containing mild and 
severe cases. We conduct numerical studies which show that, given data on sufficiently 
^ . many households, it is generally possible to discriminate between the two models by 

H I comparing the Kullback-Leibler divergence for the two fitted models to these data. 

Keywords: Household epidemic model, infector dependent severity, Kullback-Leibler diver- 
gence, multitype epidemic, varying response, final outcome data. 
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1 Introduction 

The present paper concerns models for infectious diseases in a community of households, in 
which the response to disease varies between individuals; we restrict our attention to having 
two different responses, denoted mild and severe. One common such situation is for example 
where there are asymptomatic cases showing no or hardly any symptoms but still contributing 
to the further spread of the disease. 

The reason why individuals show different symptoms may vary for different diseases. 
In the present paper we focus on two potential explanations. The first explanation is that 
the disease response is determined by individual characteristics, for example someone having 
partial immunity might become asymptomatic if infected (e.g. Staalsoe and Hviid (1998) for 
malaria and Leroy et al. (2001) in the context of ebola). The second explanation we consider 
is where the response depends on the type of infectious contact and/or whom the individual 
was contacted by. Examples where this seems to be the case are dengue fever (Mangada 
and Igarashi 1998), measles (Morley and Aaby 1997) and varicella (Mehta and Chatterjee 
2010). Ball and Becker (2006) consider the evaluation of vaccination strategies for a model in 
which infectious cases may be either mild or severe. However, their analysis is based on post- 
vaccination reproduction numbers rather than on mechanistic models such as those considered 
in this paper. 

The first explanation, where the response is determined by individual characteristics of 
the infected person, is suitably modelled using a multitype epidemic household (MT-HH) 
model (Ball and Lyne 2001). In the MT-HH model individuals are categorized into different 
types; the type of an individual may affect susceptibility to the disease and also response, 
in particular infectivity, in the event when the individual becomes infected. Quite often 
an individual's characteristics would not be known, which implies that the proportions of 
individuals of the different types in the community are unobserved. 

The second explanation can be modelled by extending the so-called infector-dependent- 
severity (IDS) epidemic model of Ball and Britton (2007) to an infector-dependent-severity 
household (IDS-HH) epidemic model. In the IDS model, the probability of an individual 
becoming a mild/severe case depends on the disease response of the person who caused that 
individual's infection. In the IDS-HH model this may also depend on whether the contact caus- 
ing the infection was a within- or between-household contact; for example, within-household 
transmission might have a higher risk of leading to severe infection. 

Once the final size distribution of the IDS-HH epidemic model is obtained, together with 
known such results for the MT-HH model, it is possible to compare the two distributions. The 
motivation for the present paper lies in this comparison. In particular we pose (and answer) 
the following question: can final size data from an epidemic outbreak with varying disease 
response be used to discriminate between the two candidate explanations for why infection 
response varies? Except in a few degenerate cases, the answer to the question is yes. This is 
illustrated numerically by showing that, in the limit as the population size tends to infinity 



in an appropriate manner, a possible outcome for either of the model is (usually) inconsistent 
with the other model. This is done by generating "data" from the MT-HH (IDS-HH) model 
and showing that the Kullback-Leibler divergence of the estimated outbreak probabilities from 
the "data" is much smaller when inference is based on the MT-HH (IDS-HH) model than when 
it is based on the IDS-HH (MT-HH) model. The final size outcome probabilities for the IDS- 
HH epidemic are obtained by numerically solving a set of differential equations, and the final 
size outcome probabilities for the MT-HH model are obtained numerically by solving a set of 
balance equations. Consequently, we have no analytical results "proving" that the two models 
are inconsistent - our arguments are instead based on numerical studies. We also consider 
data generated from finite populations and use a simulation study to demonstrate that it is 
possible to discriminate between the models using a pseudolikelihood approach. 

The paper is organised as follows. In Section [2] we define the MT-HH model and review 
final outcome results for that model. In Section [3] we define the IDS-HH model and derive an 
appropriate determinstic approximation to it. In Section H] we compare and contrast the final 
size outcomes of the two models via simulation studies, which (i) confirm the applicability 
of the asymptotic results to finite populations, (ii) strongly suggest that, as proved by Ball 
and Lyne 2010 for the MT-HH model, the final outcome of the IDS-HH model satisfies a 
central limit theorem and (iii) shed light on some interesting differences between the models. 
In Section E] we show numerically that inference from final outbreak data makes it possible 
to distinguish between the two models, using both infinite populations, as described in the 
previous paragraph, and also finite populations, where a pseudolikelihood approach (cf. Ball 
and Lyne 2010) is applied to simulated data. The paper ends with a short discussion in 
Section [HI 



2 The multitype household model 

2.1 Model definition 

The multitype household epidemic was first analysed in depth by Ball and Lyne (2001), see 
also Becker and Hall (1996) and Britton and Becker (2000). We now describe this model using 
slightly different notation. With the present application in mind, we restrict the model to two 
types and exponentially distributed infectious periods. 

Each individual is categorized as being a mild or a severe type, with the interpretation that 
if infected, the individual will become this type of infective. Additional to this, individuals 
reside in households. Let ruk^s denote the number of households having k mild and s severe 
individuals, let m„ = J2^^o ^^k,n-k denote the number of households of size n, and let m = 
^^i^n(= ^fcs'^fc.s) denote the total number of households. Further, let N = ^^i ^'^n 
denote the total population size, which is assumed to be finite. Our analysis is of the limiting 
situation in which the total number of households m tends to infinity in such a way that 
nin/m -^ pn {n = 1,2, . . .), where ^^iPn = 1 and the limiting mean household size pn = 



Yl'^=i ''^Pn is finite. It would rarely be the case that the type of an individual is known, so we 
assume that each individual is a mild case with probability Pm (and severe with probability 
1 — /3m), with the types of different individuals being mutually independent. It then follows 
that the number of mild cases in a household of size n is binomially distributed. Hence, in a 
large community, we have that TUk^s/rrik+s ~ ( fc*)/^^!! — PmY, and this holds with equality 
in the limiting situation described above. 

The disease spreads according to the following rules. Initially, a small given number of in- 
dividuals are infected (from some external force) and the remaining individuals are susceptible. 
During his/her infectious period a mild infectious individual has (global) infectious contacts 
with any given other mild individual at rate X^m/^ ^^^ with any given severe individual at 
rate \\fg/N. Similarly, an infectious severe individual has (global) infectious contacts with 
any given mild individual at rate \sm/N and with any given other severe individual at rate 
\ss f-^- Additionally, an infectious mild individual has (local) infectious contact with any 
given other mild member of his/her household at rate X^m ^^^ with any given severe mam- 
ber of his/her household at rate X^s- '^^^ corresponding rates for local infectious contacts 
of an infectious severe individual are Xgjlj and A^^. An 'infectious contact' is defined as a 
contact which results in infection if the other individual is susceptible - otherwise the con- 
tact has no effect. The infectious period of all individuals follow exponential distributions, 
with rates 7m and 75 for mild and severe infectives, respectively. All contact processes are 
governed by homogeneous Poisson processes, having rates as above. Further, all infectious 
periods and all contacts processes (whether or not either or both of the indivdiduals involved 
are the same) are assumed to be mutually independent. We assume that infected individuals 
are able to infect other individuals as soon as they have become infected, i.e. there is no latent 
period. Once an inidividual's infectious period is over, he/she recovers and becomes immune 
to further infection. The absence of a latent period, though unrealistic for most, if not all, 
human diseases, has no consequence for our present purpose, since the distribution of the 
final outcome of the MT-HH model is not changed if an almost surely finite latent period is 
incorporated (provided the rest of the model is the same). 

The MT-HH model has the following 11 parameters: 6'*^^^^ = (A]^{^, Xj^g, A^^, A^^ , X^^]^, 
-^MS' ^SM^ Xss^lM, Is, Pm)- Later we consider final size data for this model. In that situation 
we can, and hence do, assume without loss of generality that 7m = Is = 1- (The final outcome 
of a closed-population stochastic SIR epidemic of this type can be obtained by considering 
a random directed graph whose vertices are the individuals in the population and, for any 
vertices i ^ j, there is a directed edge from i to j if and only if individual i, if infected, 
has infectious contact with individual j; see, for example, Pellis et al. (2008). The set of 
people who are ultimately infected by the epidemic is given by the those individuals for 
which there is a chain of directed edges leading to them from an initial infective. Thus if, 
for example, 7m 7^ 1, we can divide all infection rates from mild infectives by 7m and then 
set 7m = 1 without changing the probability measure of the above random directed graph, 
and hence without changing the final outcome distribution. Note that this directed random 
graph explains also the above comment concerning a latent period.) It is also shown in Ball 



et al. (2004) that the 4 global infection parameters are not uniquely identifiable from final 
size data - what is identifiable are two separate linear combinations of these four parameters 
(details are given at the end of Section [2l2l) . In conclusion we hence have 7 parameters that 
are identifiable from final size data for the MT-HH model. 



2.2 Large population properties of the MT-HH model 

The MT-HH model is closely related to the model analysed in Ball and Lyne (2001). (The 
latter model allows for arbitrarily many types, non-random allocation of types of individuals to 
households and arbitrary but specified infectious period distributions.) Using essentially the 
same argument as in Ball and Lyne (2001), the MT-HH model possesses a threshold parameter 
R^,, a reproduction number for the proliferation of infected households, which determines 
whether or not an epidemic started with few initial infectives can become established in a 
large population. We now consider the final outcome of such an epidemic that becomes 
established, so implicitly we assume that i?* is above its threshold value of 1. For n = 1,2,... 
and tmiTs = 0,1,... such that tm + ^5 < n, let pn {fMi^s\d^^'^^) denote the limiting 
fraction of households of size n that have tm mild cases and rg severe cases at the end of 
an epidemic that becomes established, where the limit is as the total number of households 
m — )■ 00. An outline derivation of a method for determining p„ {TMT'rs\0'^^''^'^) is given 
below. It is a slight adaptation of the argument used in Ball and Lyne (2001), which should 
be consulted for further details. 

It is fruitful to consider first the following two-type single-household epidemic model pro- 
posed by Addy et al. (1991). Suppose that the household is of size n, that it contains k mild 
individuals and n — k severe individuals, and that all n individuals are initially susceptible. 
During the course of the epidemic, individuals avoid infection from outside of the household 
independently, with probabilities ttm and vr^ for mild and severe individuals, respectively. The 
local spread within the household is governed by the same disease dynamics as in the MT-HH 
model. Write 
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TV = (tTMjTTs) 



and denote this single-household epidemic model by E^'^'^\AS^\'k). (Recall that we as- 
sume that 7m = 75" = 1-) Let Z}^^ and Zg' denote respectively the numbers of mild 
and severe removed cases in the household at the end of the single- household epidemic, let 
p("-^)(z,j|A(^),7r) = P(4f ) = i,Zp'^ =j)iO<t<k,0<j<n-k), /^if \AW,7r) = 
E[zi"'''^] and /iJ''^'^(A(^),7r) = E[Zp''\ The probabilities p("''=)(i, j|A(^), tt) may be deter- 
mined using the following triangular system of linear equations (see Addy et al. (1991, Equa- 



tion (4))): 

(0 <ii < A;, < ji <n- A;), 
1 



where 



and 
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1 + (A; - ^l)A\^j^^ + (n - A; - ji)\l,/s 
1 



l + (A;-2i)Alf,^ + (n-A;-ji)Ag' 

The means /i]^' (A'-'^^tt) and /i^' (A*-^',7r) are easily computed once the probabihties 
p("''^)(?, j|A*^-^\ tt) have been obtained. 

Returning to the MT-HH model, suppose that there are few initial infectives, and let zm 
and zs denote respectively the proportions of individuals that are ultimately mild removed 
and severe removed, respectively. Then, if the total population size A^ is large, the probability 
that a given mild susceptible avoids global infection throughout the course of the epidemic is 

approximately ttm = exp[-(A^ZM^^ + Nzs-jf-)] = sxpI-^zm^mm + ^s>^sm)]- The corre- 
sponding probability for a given severe susceptible is vrs = exp[—{zM^MS + -^^'-^ss )]• ^^ ^^^ 
limit as the number of households m — )■ oo, the approximate probabilities tim and lis become 
exact and whether or not distinct individuals avoid global infection become independent. It 
follows that, in the event of an epidemic becoming established, the final outcome within a 
typical household of size n, that initially contained k mild and n — k severe susceptibles is 
distributed according to the final outcome of the single-household epidemic E^^'^'{A^^\7r), 
with TV = (ttm, tts) given by 

TTM = exp[-(2;MASJl^ + zsXfjlj)] and ns = exp[-{zMXMl + ^s><fj)]- (2.1) 

Note that the expected final number of mild removal cases in a household chosen uniformly 
at random is given by ^hZm- Thus, by conditioning on first the size of and then the number 
of mild individuals in such a randomly chosen household, we have that 



I^hZm 
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E/^«EurM(i-/3Mr-v£''^(A("\^). (2.2) 

n=l fc=0 ^ ^ 



A similar argument shows that 
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^^HZs = E/^«E UrM(l - /3M)"-Vf''HA^'U). (2.3) 
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n=l fc=0 



After substituting for tt from fl2.ip . equations f l2.2p and f l2.3p give a pair of non-linear equations 
for {zm, zs). These equations always have the solution {zm, zs) = (0, 0). \i R^ <1 this is the 
only solution, whilst if i?* > 1 there is (subject to very mild conditions on the parameters) a 
unique second solution in [0, 1]^, {z%j^z*g) say, giving the proportions of individuals that are 
ultimately mild and severe removed in the event of an epidemic that becomes established. It 
follows that, if tt* = (vr^,vrj) is obtained by substituting {zm^zs) = {zl^,Zg) in fl2.ip . then, 
for n = 1, 2, . . . and < vm + rs < n, 



n—rs 



pL^'^HrMyrs 



l^(MT)) ^ J2 ll]pU^-/3Mr-¥''''\rM,rs\A^'\n*) (2.4) 

k=rM ^ ' 



Calculating these final size probabilities numerically is relatively straightforward and follows 
exactly this procedure. Having substituted for tt from 02.10 . we first solve (numerically) the 
balance equations 02.20 and 02. 3p to find (2:^,2;^), substitute this into 02. ip to find (vr^,f,7rj), 
then use 02. 4p to calculate the final size distributions {pn \'CM-,'^'s\Q^^^'^'^y\- 



(G) ,(G) ,(G) x(G)n _.• r,.;„„ .* \iG) , .xiG) 



M 



Note from (^^ that any (A^^, A\^^, A^^, A^) satisfying z*mXI/^^j + zl.X^gj^ = -logvr 
and z1jX]y.jg + z^XgJ = — logvrj yields the same final size probabilities {pn {tmi ^ s\d^^^'^'')} 1 
so only these two linear combinations of (A]^^|^, X\,jg, ^smi ^ss) ^"^^ ^°^ ^^^ individual global 
infection rates are identifiable from final size data. Thus when fitting the MT-HH model we 
estimate (ttm, t^s) rather than (AJ^]^, A}^^, A^^, A^^ ). 



3 The IDS household model 

3.1 Model definition 

The infector-dependent-severity household model is an epidemic model where infected in- 
dividuals may, upon infection, become either severely infected or mildly infected, and the 
probability that an infected individual becomes mildly (or severely) infected may depend on 
both the type of its infector and whether the infectious contact is local or global. Addition- 
ally, individuals reside in households and the transmission rate is typically appreciably higher 
between individuals sharing a household. The model is defined as follows. 

Assume that there are N individuals in total, and that each individual resides in a house- 
hold. Let nin denote the number of households of size n and let m = Yl'^=i ^n denote the 
total number of households. We consider the limiting situation in which the population size 
tends to infinity in the same way as described in Section 12.11 Initially there are fc]^ mild 
infectives and k^ severe infectives, with the remaining individuals assumed to be susceptible. 
(The locations of the initial infectives is discussed later.) Mild infectives recover and become 
immune at rate 'Jm and severe infectives recover and become immune at rate 75. Thus, the 
infectious periods of infectives are assumed to follow exponential random variables, with pa- 
rameter depending on whether an infective is a mild or a severe case. While infectious, a 



(C) 

mild infective makes global infectious contacts with any given individual at rate \\j /N. If 
a contacted person is susceptible he/she becomes mildly infected with probability p^jvi ^"^^ 

(G) 

severely infected with probability 1 — Pmm'^ if ^ contacted person is already infected then 
the contact has no effect. Additionally, a mild infective has contact with any (other) given 
household member (local contact) at rate A\^ , and such a contacted individual, if susceptible, 
becomes a mild infective with probability Pmm ^^(^ severe infective with probability 1 —Pmm- 
Severe infectives have contacts according to the same rules, although with parameters A^ /N, 

PsMi -^5 ^"^^ PsM- ^11 contact processes and infectious periods are assumed to be mutually 
independent. The epidemic continues until there is no (mild or severe) infective present, when 
the epidemic stops. 

The parameters of the IDS-HH model are 9'^^^^^ = {am, -^5 ; -^m ; -^5 iPmmiPsmiPmmi 
Psm)'1m-,1s)- Note that rescaling time does not change the final outcome of an epidemic, so, 
without loss of generality, we may assume that e.g. 7m = 1, whence there are 9 parameters 
that are, in principle, identifiable from final outcome data. Note also that the directed random 
graph argument used for the final outcome of the MT-HH model fails to hold for the IDS-HH 
model, since the distribution of edges emanating from any given individual depends on the 
type of that individual, which is not determined at the outset of the epidemic and indeed 
depends on the temporal behaviour of the epidemic. Thus, fixing 7m as well as 75 would 
involve a loss of generality and the distribution of the final outcome of the IDS-HH model is 
generally not invariant to a latent period. 



3.2 Large population properties of the IDS-HH model 

Suppose that the epidemic starts at time t = and for t > 0, let X^™^,- ^ ^(t) denote the number 
of households of size n that at time t have i mild infectives, j severe infectives, k mild removed 
individuals and d, severe removed individuals. Assume now that there is a maximal household 
size nmax, so p„ = for all n > rimax- For t > 0, let X^"^\t) be the vector obtained by letting 
n, i, j, k, i vary over all possible feasible values; viz. n = 1,2,..., nmax, 0<i+j + k + i<n. 
Then {X^™'(t) : t > 0} is a density-dependent Markov population process which can be 
analysed using theory developed in Ethier and Kurtz (1986, Chapter 11). 

Suppose that m^^X^™'(0) — )> x{0) as m — )> 00, where x{0) satisfies Yln=T ^ke^'ri:o,o,k,ei^) 
< 1, so a strictly positive fraction of the population is initially infected in the limit as m — )■ 00. 
Then the above-mentioned theory of Ethier and Kurtz (1986) shows that the IDS-HH epi- 
demic process scaled by m, X{t) := X{t)/m, converges in probability to a vector of deter- 
ministic functions defined by a set of differential equations. More precisely, the component 
Xn:i,j,k/{t) = Xn:i,j,k/{f)/m = PnXn:i,j,kA't) /'^n couverges to PnXn:i,j,k/{t) defined below. The 
interpretation of Xn:ij,k,e{t) is hence the (asymptotic) fraction of the size-n households that 
at time t have i mild infectives, j severe infectives, k mild removed individuals and £ severe 
removed individuals. Using this notation we can define the (asymptotic) fraction mildly and 



severely infected at time t, iuif) and isif) respectively, by 

iM{t) = ^ ipnXn:i,j,k/(j')/fJ'H 
■n,i,j,k,e 

is{t) = ^ jPnXn:i,j,kAt)/fl'H- 
n,i,j,k,l 

The functions Xn:i,j,k,e(t) are defined by the following set of differential equations: 

x'n:i,j,k,M = [^M^PMM'^Mit) + A^ VS/^^lt) + ^M^PMMi'^ " 1) + ^f^P^llJ 

X {n- {i-l+j + k + i)) PnXn:i^l,j,kAt) 

+ \^M PMS^M[t) + As PSS ^S[t) + ^M PmS^ + ^S PsS U " ^ 

X (n- {i+j -l + k + i)) pnXn:i,j-l,k/{t) 

+ 75(j + l)pnXn:i,j+l,k/-l{t) 

- [x^S^ZMit) + Xfhsit) + Ai?^ + Af j) {n-{t+j + k + i)) 

^ PnXn:i,j,k/\t) 

- ilMi + 'ysj)PnXn:i+l,j,k-l,e{t), (3.1) 

with initial values given by Xn:i,j,kA'^) = PnXn:i,j,kA^)- 

The differential equation (13.11] applies to all relevant {n : i,j,k,i), i.e. where each of 
the indices are non-negative and i + j + k + £ < n. Vector components 'out of bounds', 
e.g. where some index is negative, are defined to be 0, for example a;3;o,-i,i,i(t) = 0. The 
first four terms in (13.11) are for households entering the state (n : i,j,k,i), explaining why 
they have a plus sign. The first term is for a household presently in state {n : i — 1, j, k^C^) 
having a mild infection and gives the overall rate for such an event to occur. The second 
term is for (n : i.,j — 1, fc,£)-households having another severe infection, the third term is 
for (n : i + l,j,k — l,£)-households having a mild removal and the fourth term is for the 
severe removals. The remaining terms describe events that cause a household to leave the 
state [n : i,j, k,i). The fifth term is the overall rate at which susceptibles in a (n : i,j, k,i)- 
households become infected, either mildly or severely and the last term is the overall rate at 
which infectives (mild and severe) in (n : i,j, fc, £)-households are removed. 

Our goal is to obtain p„ {rM,rs\9^^^^^), the limiting fraction of households of size n 
that have tm mild and rs severe cases at the end of the epidemic. If the numbers of initial 
mild and severe infectives, k^ and kg , are held fixed as m — )■ cxd, then, for large m, the 
epidemic can become established only if the household reproduction number R^, is strictly 
larger than one. (The reproduction number R^, can be obtained by approximating the process 
of infected households by a two-type branching process, the type of an infected household 

9 



being the type of its initial case; we omit the details as i?* is not required for the present 
paper.) Ideally, we would like to be able to calculate pn {rMi^s\(^''^^^^) for an epidemic 
that becomes established under these conditions. However, if k}^' and ky" are held fixed, 
then X]n=r X]fc£^n:o,o,fc/(0) = 1 and the theory of Ethier and Kurtz (1986) cannot be applied 
directly. Thus we assume instead that a very small, but strictly positive, fraction of indi- 
viduals are initially infected and approximate pn {tmi rs\0'^^^^'')-i by solving the differential 
equations ( 13.1 p numerically up to a time when the remaining fraction of infective individuals 
is negligible. More specifically, we assume that a fraction fs = 10^^ of the population is 
initially severely infective, with these infective individuals being chosen uniformly at random, 
so 

I U otherwise. 

We stop the numerical integration at the first time t' when the proportion of the population 
that is infective, i.e. iuii') + isit'), is less than 6 = 10"''' (<^ fs)- The final size probabilities 
are then given by ph, (^m,''^s|^^^'^'^^) = 3:n;0fi^rM,rsi^')- However, note that the final size 
probabilities are essentially insensitive to the initial conditions, provided that the proportion 
of index cases is sufficiently small. 

The theory of Ethier and Kurtz (1986, Chapter 11) can also be used to show that, in the 
limit as the number of households m — t- oo, the fluctuations of the stochastic model X{t) 
about its deterministic limit x{t) (defined by Xn:ij,k,e{t) = PnS:n:i,j,k,e{t)), after being suitably 
scaled, converge to a zero-mean Gaussian process, whose covariance function can, in principle, 
be determined. As in Ball and Britton (2007, 2009), this central limit theorem can be extended 
heuristically to hold also for the end of the epidemic, the time of which tends to infinity as 
r/i — )■ oo, by making a random time scale transformation in which the clock runs at rate 
m(A]^^ /a/ (t) + A^ Is{t))~^, where /m(^) and /s(t) are respectively the total number of mild 
and severe infectives present at time t in the untransformed process, cf. Ethier and Kurtz (1986, 
pp. 466-467). This yields a multivariate central limit theorem for the quantities Zn{rM,'r's) 
{n > 1, tm, rs > 0, rM+rs < n), where Z„(rM, rs) is the number of households of size n which 
ultimately have tm mild removed and rs severe removed individuals. In principle, it is possible 
to compute the covariance matrix of the limiting normal distribution numerically, though in 
practice the required computations are prohibitive except for populations comprising only 
very small households. For a population with households of sizes 1,2,..., nmax, determining 
the deterministic limit x{t) requires solving a system of n^*" = (^"™='^''+ ) _7^^^^ _ x differential 
equations and determining the above covariance matrix requires solving a system of ("^ ^~^ ) 
differential equations. For nmax = l, 2, 3, 4, 5, n^^"^ = 4, 18, 52, 121, 246, so while it is perfectly 
feasible to solve for x{t) numerically, that may not be the case for the covariance matrix. 
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4 Numerical illustrations of model behaviour 

To illustrate the asymptotic results given in the previous sections and explore some of the 
properties of the two models we have presented, we performed simulation studies of both 
models and compared some of their final size properties. In order to do this we first needed 
to select values for the parameters of our models. 

First we address these parameter choices in the MT-HH model. As mentioned in Sec- 
tion 12.11 the removal rates can without loss of generality be set to unity: 7m = 75 = 1- 
The fraction of mild types in the community was set to Pm = 0.4. This value was chosen so 
that approximately one third of all infected are mild cases (reported by Carrat et al. (2008) 
to be the case for asymptomatic cases regarding influenza). The global contact rates were 
chosen as aj^jvi — 0.25, \\,ig = 0.8, A^.^^ = 0.8 and a^jvi — ^■'^i ^^ severe infectives are more 
infectious and also mild infectives rarely globally infect mild susceptibles. The corresponding 
local contact rates were chosen as a^m — 0-2, \ms ~ ^■'^■> ^sm — 0-4 and X^m ~ ^■^' ^° 
severe infectives are also more infectious locally and mild infectives have a relatively higher 
probability of infecting mild susceptibles when compared with global contacts. The absolute 
values of the two contact matrices were chosen so that approximately 50% of the population 
becomes infected, this being a realistic value for influenza (see Ferguson et al. (2005)). The 
relative magnitude of the global and local contact rates was chosen so that both types of 
contact play a significant role in the spread of infection. 

The parameters of the IDS-HH model were chosen to be A}^^ = 1, A^. = 2, p^i/ ~ 0.8, 
pf^j = 0.2, X'-^^ = 0.5, Xf^ = 1, p55'i/ = 0.5, pf^j = 0.1, 7m = 1 and 7s = 2. These parameter 
values were chosen for the same reasons as the parameters for the MT-HH model, with the 
addition that the infectious period was set to be shorter (on average) for severe cases than for 
mild cases, having in mind asymptomatic individuals who are less likely to 'self-quarrantine' 
as they are unaware of their infection. 

The parameter common to both models is the distribution of household sizes. In this paper 
we consider two different population structures. The first is for the case where households of 
size 1, 2 and 3 are equally likely and no larger households exist, i.e. pi = P2 = Ps = 1/3, and is 
chosen largely for computational convenience. The second population structure corresponds 
to the household structure of UK in 2003 (found by typing 'household sizes' into the search 



box at http://www.statistics.gov.uk/census2001/census2001.asp ), with the simplification that 



households of size 5 and larger were truncated and all assumed to have size 5 (only 2% of 
the households had larger household size than 5, so this truncation should have negligible 
effect). The household structure for this case is given by pi = 0.29, p2 = 0.35, ps = 0.15, p4 = 
0.14, p5 = 0.07. For future reference we denote these distributions by p^^^ = (1, 1, l)/3 and 
p(5) — (29, 35, 15, 14, 7)/100. We also note that 3 is the minimum value of nmax for both 
models to be, in principle, identifiable. 

For both models we ran 10,000 simulations of systems with 10,000 households, the sizes 
being given by p = p^^\ (We treated p as giving the proportions of households of different sizes, 
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rather than having random household sizes with distribution p.) In order that minor outbreaks 
are unhkely, we initiated the epidemics with 10 infectives, each randomly chosen in different 
households of size 5; in the MT-HH model these individuals may be mild-type or severe-type 
(with respective probabilities /3m and 1 — (3m) and in the IDS-HH model we specified that 
they are all severe cases. For simulations that result in more than 0.15 of the population 
becoming infected we then recorded the overall final size amongst initial susceptibles and 
the within-household final sizes amongst households that are initially completely susceptible. 
(Inspection of histograms (not shown) of the final proportion of individuals infected suggests 
that this cutoff is appropriate for separating minor and major outbreaks.) 

Figures [1] and [2] show histograms of the numbers of individuals ultimately mildly and 
severely infected in, respectively, the 9,992 simulations of the MT-HH model and the 9,993 
simulations of the IDS-HH model that resulted in major outbreaks. Overlaid on these his- 
tograms are probability density functions (scaled so as the area under them matches that 
of the histograms) of normal distributions with the same mean and variance. The excellent 
agreement between the histograms and density functions in Figure [1] is expected in view of 
the central limit theorem for the MT-HH model of Ball and Lyne (2001) and in Figure H] 
this lends credence to the central limit theorem discussed in Section 13.21 above for the final 



outcome of the IDS-HH model. Though the mean values of these distributions are similar 
(/iJvf ^^ ~ 4,525, pSi^^^ ^ 4,854; plf ^^ ^ 9,835, /il.^^'^^ ^ 10,008) - indeed the parameter 
values were chosen with this intention - it is interesting to note that the variability is rather 
different in the two models. The spread of the distribution of the number of mild cases in 
the IDS-HH model is appreciably larger than that in the MT-HH model {a\j ^ 93 and 
aj^j ~ 150) and, though not to the same extent, the distribution of the number of severe 
cases is also slightly more spread in the IDS-HH model {a^ ~ 167 and a^ ~ 218). Part 
of the reason for this is that in the MT-HH model the types of individuals are determined 
in advance, but in the IDS-HH model the types of the infected individuals depend on the 
evolution of the epidemic and some feedback may occur (though with different parameters it 
might potentially be positive or negative). 

Tables [1] and [2] give further information about the within-household outcomes of major 
outbreaks in the two models. These tables give, for each household size n, estimates from the 
simulations of the probability that a typical individual in a household of size n is ultimately 
(i) mildly infected (pm), (h) severely infected (ps) and (iii) infected (pinf = Pm +Ps), and 
also of the probability that a case in a household of size n is severe {ps/pinf)- The figures 
in parentheses are the corresponding infinite population asymptotic quantities obtained from 
{pr''\rM,rs\e^^''^'^)} and {p'n'''\rM,rs\e('''''^)}, respectively. 

In both cases we observe good agreement between the deterministic and estimated stochas- 
tic quantities. Also observe that, in both models, the proportion of individuals infected in- 
creases with household size n. This a consequence of local spread being greater in larger 
households. In the MT-HH model the proportion of cases that are severe decreases with n, 
whereas this proportion increases with n in the IDS-HH model. In the IDS-HH model this 
simply reflects the fact that local infections are very likely to result in severe cases. In the MT- 
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Figure 1: Histograms of final outcome of major outbreaks in simulations of the MT-HH model 
in a community of 10,000 households, with matched normal approximations superimposed. 
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Figure 2: Histograms of final outcome of major outbreaks in simulations of the IDS-HH model 
in a community of 10,000 households, with matched normal approximations superimposed. 



HH model, however, an individual's type is determined in advance rather than by the spread 
of infection and this results in a 'saturation effect' of sorts. Note that the proprtion of cases in 
households of size 1 that are severe (0.7191) is larger than the proportion of individuals that 
are of severe type (0.6). In larger households more local spread is expected than in smaller 



MMi -^MSi 



,W 



A^^^ and X^g are independent of household size) and the rates are 



households (as A 

such that severe types are more likely to be infected. Indeed, in a very large household we 
would expect everyone to be infected, in which case the proportion of cases that are severe 
must be equal to the proportion of individuals of severe type. For any household size, the 
proportion of globally contacted individuals that are severe is 0.7191. However, since local 
spread increases with households size and there is greater scope for local infection amongst 
mild types than severe types (as relatively fewer are infected globally) the proportion of cases 
that are severe must decrease with household size. 
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Table 1: Properties of MT-HH epidemics that become established. 



n 


Pm 


Ps 


PiNF 


Ps/PmF 


1 


0.1272 (0.1273) 


0.3255 (0.3256) 


0.4527 (0.4529) 


0.7190 (0.7189) 


2 


0.1585 (0.1585) 


0.3751 (0.3753) 


0.5336 (0.5337) 


0.7030 (0.7031) 


3 


0.1923 (0.1925) 


0.4228 (0.4229) 


0.6151 (0.6154) 


0.6873 (0.6872) 


4 


0.2270 (0.2271) 


0.4656 (0.4658) 


0.6926 (0.6929) 


0.6723 (0.6722) 


5 


0.2602 (0.2603) 


0.5020 (0.5021) 


0.7621 (0.7624) 


0.6587 (0.6586) 



Table 2: Properties of IDS-HH epidemics that become established. 



n 


Pm 


Ps 


PlNF 


Ps/PlNF 


1 


0.1815 (0.1822) 


0.2870 (0.2865) 


0.4685 (0.4687) 


0.6126 (0.6113) 


2 


0.1969 (0.1976) 


0.3546 (0.3542) 


0.5515 (0.5517) 


0.6430 (0.6419) 


3 


0.2095 (0.2104) 


0.4267 (0.4261) 


0.6362 (0.6364) 


0.6707 (0.6695) 


4 


0.2190 (0.2196) 


0.4980 (0.4975) 


0.7169 (0.7171) 


0.6946 (0.6937) 


5 


0.2229 (0.2250) 


0.5671 (0.5638) 


0.7901 (0.7888) 


0.7178 (0.7147) 



Clearly the above phenomena depend on the parameter values chosen in the two models. 
For example, in either model, simply interchanging the labels of the two types results in the 
opposite effect of household size on Ps/pinf being observed. 



5 Model discrimination 

Suppose data from a population with household structure given by {pn} are generated from one 
of the two models, with some given parameters ^(^^) or 9^^^^\ as appropriate. An important 
inference, or discrimination, problem in light of two possible models is then whether it is 
possible to determine which of the models the data come from. It is hard to give an analytical 
answer to this question since the final size probabilities are not explicit. We address the 
question with a numerical investigation. For our purposes, the data are the distributions of 
within- household final sizes q = {qn{rM,rs), < vm + rs < n, < n < nmax}- We consider 
the case where q is the asymptotic (m — )■ oo) final size distribution derived from one or other of 
the models using the methods described in Sections 12.21 and [ 3.21 in order to determine whether 
the two models actually produce different final size distributions. We also consider the case 
where q is derived from stochastic simulations of one or other of the two models (i.e. with 
m finite), to determine whether or not any difference between the two models is sufficiently 
pronounced to be detectable with a dataset that resembles more closely one available in real 
life. 
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In the remainder of this section we describe first, in Section [3?T| how we generate the data 
that we use, both from infinite and finite populations, then discuss, in Section [3^ how we fit 
the models to a given final size distribution. In Section 15. 3[ we describe our main findings 
concerning whether the MT-HH and IDS-HH models can be distinguished on the basis of final 
size data. In Section 15.41 we motivate our use of the KuUback-Leibler divergence as a tool for 
model fitting and discrimination and finally, in Section 15. 5[ we discuss identifiability issues 
that arise in fitting the models to final outcome data. 



5.1 Data generation 

Final size data for an infinite population are generated using the methods described in the 
previous sections. For the MT-HH model we solve equations (2.1)-(2.4) numerically and 
for the IDS-HH model we solve the differential equations (13. ip numerically, as described in 
Sections 12. 21 and [3. 21 respectively. To generate final size data for finite populations we simulate 
an outcome of the relevant stochastic process according to the model description in Section [2?T] 
or 13.11 as appropriate, (with 10 initial severe infectives, each in separate households of size 
'^max; to iucrcasc the chance of a major outbreak occuring). If a major outbreak does occur 
(which we take to be more than 0.15 of the population becoming infected) then we calculate 
the empirical final size distribution considering only the households in that simulation that 
had no initial infectives. In either case we denote by g = {qn{i^M,rs)} the 'target' household 
final size distributions that we try to reproduce from the model we choose to fit to the data. 



5.2 Model fitting 

We now describe the algorithm we use to fit each model to given final size data q = {gn(^M, f^s)}- 
The goal is to find parameters 6 of the model we are fitting so that the distance between the 
final size distributions {pn{'i"Ai,rs\9)}, corresponding to 9^'^'^' or 6^^^^\ and the 'target' final 
size distributions q = {qn{rM, rs)} is as small as possible. We measure this distance using the 
KuUback-Leibler (K-L) divergence 

m = DKL{qM9)) = Y.Pr^ll gnK/, r5) log Pp^ . (5.1) 

The use of the K-L divergence is motivated by its well-known relationship with likelihood- 
based inferential procedures (see, for example. Bishop et al. (1975, pp. 344-348)), which is 
discussed in more detail in Section [531 We minimise f{6) numerically using Matlab's fmincon 
constrained optimisation routine. This requires selecting a starting point ^o for fhe parame- 
ters; we choose these starting values independently at random, the rate parameters from an 
exponential distribution with mean 1 and proportion/probability parameters uniformly from 
the interval (0, 1). In the case of fitting the IDS-HH model we find that the numerical opti- 
misation is more difficult and that it is beneficial to sample several (we use 20) such possible 
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starting points 6q and then start the numerical optimisation routine at the best of these points 
(i.e. that with smallest /(6'o)), so that the numerical routine is more likely to start at a point in 
parameter space that is at least moderately compatible with the target final size distributions 
q. We describe this process of choosing a starting point for and then running the optimisation 
routine as a single 'run' of our algorithm (i.e. model fitting procedure). 

Calculating p{9) = {Pn{^M,rs\(^)} using the methods described in Section fI72\ or IX^ is 
straightforward and in principle evaluating f{6) is then trivial as long as p has no zero entries, 
i.e. as long as the parameter vector 6 results in the model being super-critical. However, there 
are numerical problems that can arise when calculating the K-L divergence as in equation (15. ip . 
These problems arise due to so-called 'catastrophic cancellation' which occurs when using the 
formula (15. ip if q and p differ only slightly. The terms g^ log(gj/pj) are all small (since p and 
q are close) but are of differing signs (since sometimes g^ > Pi and sometimes vice- versa), 
thus when the sum ^^ qi log{qi/pi) is close to zero there can be catastrophic cancellation 
and the calculated sum can be wildly inaccurate. We resolve this by using the Taylor series 
approximation slog(s/t) ^ (s — t)^/2t about s = t (cf. Bishop et al. (1975, Lemma 14.9-1)), 
which implies 

f{0) = DKLiqMOj) - y^Pn y 7^—, ^T • 5.2 

n=l rM,r.s 

This approximation becomes exact as p — ?■ g so using it when the calculated K-L distance 
is small gives a good approximation and avoids numerical problems. Numerical experiments 
comparing the calculated values of f{6) using (15. ip and (15. 2p show good agreement, improving 
as f{9) becomes smaller (precisely as expected), but when f{9) is less than about 10~^ we 
begin to see significant disagreement. Therefore, all our calculations of K-L distance initially 
use (15.11) but if the result is smaller than 10~^ we recalculate using (15.21) . 

The random starting values of 9^^'^^'' and 9^^^^^ in our fitting procedure will sometimes 
be poor (i.e. give large values of f{9)) and result in the optimisation routine staying in a part 
of parameter space that gives a very poor fit. Thus, when fitting a model to data we run our 
algorithm many times over to ensure that as much as possible of the parameter spaces are 
explored. The number of these runs necessary is somewhat variable; this issue is addressed in 
Section 15.51 

Initially we focus simply on the smallest of the K-L distances f{9) of the model from the 
data that we find for each combination of dataset and model. In Section 15.51 we explore in 
more detail the variability of the f{9) from run to run of our algorithm and also examine the 
behaviour of the corresponding parameter estimates 9. 
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Table 3: Best fits of eacli model to final size distributions obtained from all four combinations 
of household size distribution and model. 



.(3) 



data 



P = P' 
MT-HH IDS-HH 



3.4x10-11 1.5x10-3 
4.7x10-5 8.9x10-9 



P 
MT-HH 



P 



(5) 

IDS-HH 



2.0x10-11 6.8x10-3 
1.1x10-^ 3.3x10-* 



1^ MT-HH 
§ IDS-HH 



5.3 Model discrimination 

5.3.1 Infinite data 

To determine whether or not each model is capable of producing the final size distributions 
generated by the other model we fit a given final size distribution to both models, in the 
expectation that the correct model will fit appreciably better. We find that the correct model 
can be made to fit as well as we please by tightening the stopping criteria of the numerical 
optimisation routine but that there is a definite non-zero lower bound for f{6) when we fit 
the wrong model. Further details of this are given in Section 15.51 

We summarise our findings by way of Table [3l which shows the best fits obtained from 
100 runs of our algorithm (measured by f{6)) obtained when fitting the IDS-HH and MT-HH 
models to the (asymptotic) final size distributions produced from each of the models (with 
parameter values as in Section Hj) with each of the household size distributions p*^^) ^^^1 p'-^^ . 
Table [3] demonstrates the significant differences in fit obtained when fitting the two models to 
each data set (i.e. each column of the table). In Tableland the following discussion, 'data' 
refers to the model that generated the given final size distributions we fit to and 'model' refers 
to the model we fit to these data. 

It can be seen that using the household size distribution p^^'' which includes households 
up to size 5 makes no qualitative difference to these conclusions. However, it is interesting 
to examine the effect of larger households on the above K-L distances. Table H] shows the 
contribution to the best final K-L distances in Table [3] from households of each size, which 
amounts to separating out the summands pn ^rM,rs ^n{rM, rs) log(g„(rAf , rs)/pn{rM, rs\0)) in 
equation (15.11) . The breakdown of the best final K-L distances suggests that an apprecia- 
bly greater contribution to the K-L distances comes from larger households than would be 
expected simply based upon the proportions of households of different sizes present. This 
perhaps suggests that data collection effort might be focused somewhat more on larger house- 
holds; though of course this depends crucially on our assumption that the same transmission 
parameters apply in households of all sizes. 
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Table 4: Breakdown of contribution to final K-L distances by households of different sizes 
when p = p^^\ 





MT-HH model, 


IDS-HH model, 


MT-HH model. 


IDS-HH model. 


n 


MT-HH data 


IDS-HH data 


IDS-HH data 


MT-HH data 


1 


4.4x10-12 


7.6 X 10-1° 


2.0x10-5 


2.0 X 10-^ 


2 


1.1 X 10-11 


1.0 X 10-9 


3.2 X 10-5 


1.1 X 10-5 


3 


1.8x10-11 


7.1 X 10-9 


1.4x10-=^ 


3.6 X 10-5 


total 


3.4x10-11 


8.9 X 10-9 


1.5 X 10-3 


4.7 X 10-5 



We see (Table E]) that the final size distribution generated by each model using somewhat 
realistic parameter values cannot be captured by the other model. To investigate whether 
this conclusion holds for the models in general, we need to do this comparison for a range 
of (supercritical) parameter values. We expect that the fits will be poor except possibly for 
some degenerate cases where the models can produce the same final size distributions. 

To test whether one model can reproduce final size data from the other, we select param- 
eters for one model at random, resampling if the corresponding R^ < 1, and calculate the 
corresponding final size distribution, then fit the other model to this 'data'. When selecting 
the random model parameters to use, each parameter is chosen independently, rate parame- 
ters from an exponential distribution with mean 1 and probability parameters uniformly on 
[0, 1]. We then repeat this experiment many times so that we explore a range of parameter 
combinations of the model from which we derive our data. The final K-L distance (best fit) 
f{9) that we report for each paremeter combination is the best fit obtained in 5 runs of our 
algorithm. (When fitting to final size data that is not exactly reproducible by the model we 
are fitting we find that the variability of f{9) between runs is very small and thus 5 runs is 
more than sufficient to be confident that we have found the best-fitting model; see Section 15.51 
for further details.) 

Figure [3] shows histograms of the best K-L distances f{6) obtained when fitting one 
model to final size data generated by the other with random parameters and household size 
distribution p^^) We have fitted the MT-HH model to 10,000 random IDS-HH datasets but 
owing to the computational expense of fitting the IDS-HH model we have only fitted it to 
300 random MT-HH datasets. It can clearly be seen that for most parameter combinations 
the final size distributions cannot be reproduced by the wrong model. Of course the correct 
model can reproduce these final size distributions and to confirm this we have also fitted the 
correct model to many of these data. As expected, the correct model fits appreciably better 
than the wrong model except in the degenerate cases discussed below when both models can 
be fitted to the data (details not shown). 



(a) fitting MT-HH model to IDS-HH data 



(b) fitting IDS-HH model to MT-HH data 
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Figure 3: Smallest K-L distances f{6) obtained when fitting the one model to the final size 
distributions from the other with random (super-critical) parameters. 



Further analysis of the cases where the 'wrong' model fits the data relatively well {f{0) < 
10^^) reveals at least one of the following reasons. In either model, if the process is only just 
super-critical, i.e. R^, is only slightly larger than 1, then many of the quantities qn{rM,fs) 
in (15.11) are very small so relatively fewer of the summands contribute to the sum and it is 
somewhat easier for the wrong model to be able to fit the data. 

In the IDS-HH model, two further situations arise where the wrong (i.e. MT-HH) model 
fits the final size data quite well: (i) one or both types of individual makes very few local 
contacts, i.e. min{A 



(L) WL) 



M ' 



^^s^' /is} is close to (recall 7a/ 



1) and (ii) local contacts by mild 
and severe individuals are approximately equally likely to cause mild/severe cases, i.e. \P]^\^ — 
Pg^^l is close to 0. In case (i), within-household spread essentially only involves one type of 
individual making contacts and the local infection rate and probability can be tuned to produce 
almost any local final outcome distribution. In case (ii) local infection processes become like 
those in the MT-HH model because each (locally infected) individual becomes mild or severe 
(independently) with the same probability p^^ ~ Psm- -'■^ ^^^ MT-HH model there are also 
two further situations where the wrong (i.e. IDS-HH) model can reproduce the final outcome 
data well. The first is if there is essentially only one type of individual, i.e. (3m is close to or 
1; if one type is not present it is trivial that the two models coincide (in the sense that they 
can produce the same final size distributions). The second case is where the disease is highly 
globally infectious amongst one type of individual, i.e. tt has (at least) one element close to 0. 
Here local transmission is essentially a one-type process and again the models coincide. 



5.3.2 Finite data 

We have just seen that it is possible to discriminate between the two models using final size 
data from an infinite population. Real data of course never pertains to an infinite population, 
so in the present subsection we perform the same type of analysis except that data are now 
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(a) simulated IDS-HH data 



(b) simulated MT-HH data 



-IDS-HH model 
-MT-HH model 




-MT-HH model 
- IDS-HH model 




Figure 4: Smallest K-L distances f{9) obtained when fitting both models to the output of 25 
separate empirical final size distributions from simulations of each model. 



generated from the stochastic models in a community oi m = 10, 000 households. The data, 
generated both from the stochastic MT-HH model and the stochastic IDS-HH model, hence 
consist of empirical final size distributions rather than the exact asymptotic distributions. 
The model fitting procedure is exactly the same as before but we now use the empirical final 
size distribution as the target final size distribution q in (15.11) . 



From each model, with parameter valules as in Section H] and household size distribution 
p(3) J we generated 25 independent empirical final size distributions taken from simulations on 
systems of 10,000 households that resulted in a major outbreak and then fitted both models 
to each empirical final size distribution. Figure S] shows plots of the fit of each model to each 
dataset, the fit being measured by the smallest value of f{9) found in 5 runs of our algorithm. 
(When fitting to empirical final size distributions we find that the variability of f{6) between 
runs is very small and thus 5 runs is more than sufficient to be confident that we have found 
the best-fitting model; see Section [531 for further details.) For clarity, the results in the figure 
have been ordered according to the best fit of the true model. 

From these plots it is immediately clear that the correct model (i.e. the one that generated 
the data) has the best fit on most occasions (24 out of 25 for the IDS-HH data and 22 out of 
25 for the MT-HH data). It also seems clear that it is generally easier to rule out the MT- 
HH model when looking at data from the IDS-HH model than vice-versa (the gap between 
the two lines is generally much larger in plot (a) than in plot (b)). Especially intriguing is 
the observation that, when the data are from the MT-HH model there is a clear association 
between the K-L distances to the best IDS-HH and MT-HH model; however there is much 
less, if any, association when the data are from the IDS-HH model. This may be an artifact of 
the fact that, as can be seen from Figure [3|, the IDS-HH model is generally able to fit MT-HH 
data better than the MT-HH model can fit IDS-HH data. 
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5.4 Pseudolikelihood motivation for use of K-L divergence 

In this subsection we motivate our choice of the Kullback-Leibler divergence for assessing the 
distance between the two models by relating it to a maximum pseudolikelihood estimation 
procedure. For ease of presentation our arguments are informal, rather than fully rigorous. 
Suppose, as above, that we have data {qn{rM,^s)} from an epidemic in a community of m 
households. If we make the approximation that the outcomes in different households are 
mutually independent then the likelihood of these data under one or the other of our models 
is given by 

HO) = n n [PnirM,rs\e)r-''-^'-^"'-'^\ (5.3) 

n=l rM,rs 

where m„ is the number of households of size n in the community, and [6, Pn{fM,^s\(^)) is 
either (^(*^^), pit'^\rM,rs\e^''^^)) or (0(^^^), pli'''\rM, 

rs\9^^^^^)). In reality, (15. 3 p is a pseudolikelihood since the outcomes in distinct households are 
dependent, as they are part of the same community-wide epidemic, though the dependence is 
small (the covariance of the final outcomes in distinct households is of order l/m for large m; 
cf. Ball and Lyne (2010)). The maximum pseudolikelihood estimator of 6, denoted by 6, is 
obtained by maximising L{6), or equivalently by maximising l{6) = log L[6), which is given 
by 

l{9) =m'^pn ^ qn{rM,rs) log Pn{rM,rs\9). 

Note that maximising l{6) is equivalent to minimising the Kullback-Leibler divergence 
Dkl{q\\p{9)), defined by (15. ip . Moreover, the pseudolikelihood ratio goodness-of-fit test statis- 
tic. Am say, for assessing the adequacy of the model for these data is given by 

- 2 log Am = 2mDKL{q\\p{e)); (5.4) 

cf., for example. Bishop et al. (1975, Equation 10.2-6), who consider testing the goodness-of-fit 
of a specified multinomial model. 

Now, for example, suppose that these data {qn{.fM^fs)} were actually generated by the 
IDS-HH model with parameter 9^^^^\ but that we fit the MT-HH model. Then (cf. Sec- 
tioning]) Qni'^^Mi fs) — > Pn (fM, rs|6'*^^^'^^) as m — > oo, whence 6**^*^^) — > 6; ^ as m — )■ oo, 
where #^^^ minimises DKL{p^^^'^\e^^^^^)\\p^^^^\e^^^^^)) with respect to O^^^^l (Here, ^ 
denotes convergence in probability.) Hence, using (15. 4p . 

--21ogAm^2DxL(p('^^H^('^^))lb('''^H^r^))) as m^oo. 
m 

If, instead, we fit the IDS-HH model with parameter 9'^^^^\ then 9'^^^^^ -^ 9^^^^^ as 
m^oo and, since D^i(j9(^^^)(0(^^^))||p(^^^)(^(^^^))) = 0, -2m-' 
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log Am — > as m — )■ oo. In these circumstances, — 21ogAm asymptotically equals the usual 
chi-square goodness of fit test statistic 



= 1 i-M^rs 



m„rf^^VM,r5|^(^^^)) 



However, dependencies between the households imply that X^ may not have the usual asymp- 
totic x^ distribution; instead the asymptotic distribution of X^ is a linear combination of d 
independent Xi random variables, where d is the degrees of freedom of the usual chi-square 
test (cf. Ball and Lyne (2010)). Nevertheless, (15.41) gives a guide for interpreting both our in- 
finite and finite population model discrimination results. Moreover, if these data {qnijM, rs)} 
come from a small fraction, Sq say, of the households among which the epidemic is spreading, 
as is often the case in practice, then, if the model is correct, the asymptotic distribution of 
X"^ is very close to the usual Xd distribution, the approximation being exact in the limit as 
£o ; (cf. Ball and Lyne (2010)). 

Clearly there are precisely analagous results which hold if the data instead come from the 
MT-HH model. 



5.5 Identifiability and model fitting 

In this subsection we investigate how our model fitting methodology works in practice. We 
find that there are two key issues that infiuence the overall behaviour of our algorithm. The 
first is the striking difference in the distribution of final K-L distances f{9) in the situations 
where the target final size distribution q can or cannot be captured exactly (to numerical 
accuracy) by the asymptotic version of the model we try to fit. The target final size distribution 
cannot be captured exactly when either (i) we try to fit the wrong model or (ii) the data is 
from a finite population. (This has the important consequence that when fitting a model to 
empirical final size distributions we need only run the algorithm a few, say 5-10, times to 
be confident that we have found the best possible fit.) We therefore restrict our attention 
here to an exploration of seeking to fit the models to 'data' which are the asymptotic (m — )■ 
oo) distributions corresponding to the parameter values given earlier, with household size 
distribution p^'^^ = (1, 1, l)/3. We refer to these data with this household size distribution and 
the parameters given previously as g(^^) and q(^^^\ In the course of this we clearly see the 
second issue that arises, namely that there are identifiability issues in the IDS-HH model. In 
the IDS-HH model as we parameterise it, it seems that some parameters are identifiable while 
some are more difficult to identify, though we can find functions of these parameters which 
do appear to be identifiable. 

5.5.1 Fitting the correct model 

Firstly we look at fitting each model to data generated from that same model, so we should 
be able to recover the input parameters used to generate the data and find that f{6) is very 
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Figure 5: Profiles of final K-L distances obtained when fitting both models to both g(*'^^) and 

{IDS) 



Q 



close to 0. Figure E] shows density estimates (essentially smoothed histograms, which we use 
for ease of display) of f{6) for the best 90 of 100 runs of our algorithm when fitting each model 
to data generated by that model. We use only the best 90% of runs so as to exclude the poor 
fits sometimes obtained for the reasons explained in the second paragraph of Section 15.21 (For 
comparison, Figure [5] also shows the smallest f{6) values found when fitting each model to the 
data generated by the other model; these are displayed as points rather than densities since, 
as shown in Section [5.5.31 in these circumstances the variability of f{6) is very small.) This 
figure shows that our algorithm consistently finds model parameters 6 which quite accurately 
reproduce the target final size distributions, indicated by the very small values of f{9). We 
see shortly why the MT-HH model can be fitted to its own final size distribution rather better 
then the IDS-HH model. If we examine the parameter estimates 6 that yield these final K-L 
distances we see (Table E]) that the MT-HH model recovers the parameters used to generate 
qi^T) ^[-^i^ g^ high degree of accuracy and very little variability, whereas when we fit the 
IDS-HH model to q^^^^^i we find (Table [6]) that several parameters are estimated quite poorly. 

Table 5: Summary of parameter estimates when fitting the MT-HH model to g(^^) (best 90 
of 100 runs). 



Parameter 



ttm 



TTS 



A 



(L) 

MM 



X 



(L) 
MS 



^SM 



^SS 



/3. 



M 



True value 

Mean 
Std. dev. 



0.7263 0.5224 
0.7263 0.5224 
0.00003 0.00004 



0.2000 


0.4000 


0.4000 


0.8000 


0.4000 


0.2000 


0.4000 


0.4000 


0.8000 


0.4000 


0.00004 


0.00013 


0.00008 


0.00022 


0.00004 
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Table 6: Summary of parameter estimates when fitting tfie IDS-HH model to q^^^^^ (best 90 
of 100 runs). 



Parameter 



,{G) ,{G) ,{L) ,{L) (G) (G) (L) (L) 

^M ^5 ^A/ ^5 PmM PsM PmM PsM [S 



True value 

Mean 

Std. dev. 



1.0000 2.0000 0.5000 1.0000 0.8000 0.2000 0.5000 0.1000 1.5000 
1.7805 3.9959 0.5028 4.2510 0.3924 0.4954 0.4935 0.0932 8.5436 
0.6781 3.4214 0.0017 2.4507 0.1662 0.2507 0.0040 0.0046 4.9344 



5.5.2 Identifiability in the IDS-HH model 

The poorer recovery of input parameters in the IDS-HH model can to a large extent be 
explained by issues of identifiability. We mention in Section |2] that it is known that in the 
MT-HH model the global rate parameters are not uniquely identifiable from final size data 
(Ball et al. 2004) and for this reason we estimate the probabilities tt rather than the global rates 
{■^MMi ^MS^ ^SMi '^ss)- However, the IDS-HH model we propose is new so such identifiability 
issues have not been explored. Moreover, identifiability is difficult to study rigorously for this 
model as there is no analytical expression for the household final size distributions. Careful 
examination of the parameter estimates when fitting the IDS-HH model to q^^^^^ suggests 
that some identifiability issues are present here. In particular, in our parameterisation of the 
IDS-HH model there are three combinations of parameters that seem identifiable whilst some 
of the individual parameters are very difficult to identify separately. 

The first of these combinations is A^ and 75; our algorithm estimates the ratio A^ /"Js 
extremely well (see Table [7]), but has difficulty identifying the precise values of these param- 

(G) 

eters. The second set of troublesome parameters consists of the global contact rates X)^ and 

(G) 

Xg and the removal rate 75. If the removal rates 'Jm and 75 are known then the relationship 
ttg = exp{-{zMXfj^ /-fM + zsXf^ /'js)), where ttg = gi(0,0) (= ^g„(0,0) for any n < n^ax) is 
the probability that a given individual avoids global infection, specifies a linear equation that 

(G) (G) 

X\j and Xg must satisfy. If we assume that the removal rates are both known then our algo- 
rithm identifies the correct linear combination of global contact rates very easily but finds it 
very difficult, though possible, to find the most likely values of these parameters individually. 
However, we assume that (one of) the removal rates is unknown and, as just discussed, not es- 
timated very well; thus the global rates are generally not estimated very reliably. Nevertheless, 
when the initial guess for 7s is close to its optimum (correct) value, A^ and the above linear 
combination of X\^ and A^ are also estimated easily and reasonably accurately, and a very 
good fit is obtained. In the latter situation it is also possible to recover the individual rates 
X\j and A^ with our algorithm but this is far more difficult. That zm^m hu + ^s^s hs 
is estimated well is demonstrated in Table [71 in which zm and zs are given by their observed 
values in the (infinite) data. 
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(G) 

The other parameters with identifiabihty issues are the global infection probabilities Pmm 

(G) 

and PsM- Some information about these parameters can be obtained by considering households 
of size 1 which become infected. Focus on such a household and suppose that there are in 
total Ym mild and Ys severe infectives in the population just prior to its infection. Then the 
probability that this infection is mild is 

^M^M PmM + ^S-^s PsM 

YmX^S^+YsX^ 

Now, Ym and Ys are random and vary throughout the epidemic. A crude approximation 
is to replace the ratio Ym/Ys by ^m^zm/'Js'^zs, the latter taking into account the different 
infectious periods of mild and severe infectives. Thus the proportion of infected households 
of size 1 that are mildly infected is approximately 



zm^m PmmIim 



^s^f^pfuhs 



AG) 



>(G) 



zm^m /iM + zsXs hs 



leading to the relationship 



p1^^"^(i,o|^(^^^)) 



SG)AG) 



iG)iG) 



zmXI/pI/m/Im + zsX's PsmHs 



p1^"^^(i,o|^(^^^))+pP^^(o,i|^(^^^)) 



{IDS). 



AG) 



AG) 



zmX'-m llM + zsXg ' hs 



(5.5) 



We have seen above that the denominator in the right hand side of fl5.5p can be estimated 
well, hence it is reasonable to expect that the numerator might be too. That this is indeed 
the case is borne out in Table [71 



Table 7: Functions of estimated IDS-HH model parameters when fitting IDS-HH model to 
g(^^^) (best 90 of 100 runs). 



Function 


ZAiXij^/^fM + zsXi. Vts 


A?V7. 


ZNlXf^P^S^^/lM + ZsXf^pfijhs 


True value 

Mean 
Std. dev. 


0.50669 
0.50672 
0.00003 


0.50000 
0.49807 
0.00113 


0.21340 
0.21069 
0.00170 



5.5.3 Fitting the incorrect model 

We now turn our attention to the situation where we try to fit one of the models to final 
size data arising from the other model. Fitting the MT-HH model to IDS-HH data gives 
parameter estimates summarised in Table 13 While there is more variation in the MT-HH 



25 



parameter estimates than when we fit to data from the MT-HH model, the variation is still 
relatively small. Furthermore, the variation in the final K-L distances f{6) is very small 
(mean 1.46 x 10"^, st. dev. 1.5 x 10^^'^), giving confidence that (i) we have found the region 
of parameter space where the MT-HH model can best reproduce the data from the IDS-HH 
model and (ii) that the best fitting MT-HH model does not reproduce the IDS-HH final size 
distribution very well. For comparison the minimum of these K-L distances is also shown in 
Figure as a point rather than a density because the variation is so small. 



Table 8: Summary of parameter estimates when fitting the MT-HH model to q^^^^^ (best 90 
of 100 runs). 



Parameter 



TTAf 



T^S 



A 



MM 



A 



MS 



A 



SM 



A 



(L) 

ss 



/3. 



M 



Mean 
Std. dev. 



0.5210 0.6450 1.3712 0.2561 0.0509 0.8990 0.3373 
0.00003 0.00001 0.00035 0.00003 0.00003 0.00007 0.00002 



Lastly we consider fitting the IDS-HH model to the MT-HH data; see Tableland Figured 
Here we see variations in 6 roughly the same as those seen when fitting the IDS-HH model to 
data it can reproduce exactly. Again we find that the estimates of X\j and the p^^'^s show little 
variation and we also see the same identifiability issues present. Although the estimates of 75, 
Xg , Xg and the p'^'^^'s individually vary wildly we find that A^ /7s, ^mA'-'^Vtm + ^^sA^ 775 
and zmX)^ p\^j^/'jM + zsXg Psm/1s show very little variation (see Table [10]). Similarly to 
when we fit the MT-HH model to the IDS-HH data, we find very little variability in the final 
K-L distances that we find, for the 90 smallest values the mean and st. dev. are 4.69 x 10^^ 
and 1.0 X 10~^, respectively. Though these K-L distances certainly seem bounded away from 
zero, suggesting that the IDS-HH model cannot reproduce the MT-HH final size distribution, 
they are appreciably smaller than when fitting the MT-HH model to IDS-HH data. 



Table 9: Summary of parameter estimates when fitting the IDS-HH model to g(*^^) (best 90 
of 100 runs). 



Parameter 



A 



(G) 
M 



X 



(G) 



A 



(L) 
M 



X 



S 



(G) 

Pmm 



(G) 

Psm 



Pmm 



Psm 



Is 



Mean 
Std. dev. 



2.3496 4.0997 0.1982 4.2774 0.2351 0.5028 0.2757 0.3470 9.7983 
0.6919 2.2107 0.0001 2.8475 0.1548 0.2579 0.0038 0.0037 6.5202 
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Table 10: Functions of estimated IDS-HH model parameters when fitting IDS-HH model to 
g(*^^) (best 90 of 100 runs). 



Function 



,(G) 



(G) 



zm>^m hM + zsXg hs 



^fhs 



l(G)JG) 



iG)JG) 



zm^m'Pmm/1m + zsX^s PsmMs 



Mean 
Std. dev. 



0.50504 
0.00003 



0.57068 
0.00023 



0.13909 
0.00143 



6 Discussion 

In this paper we define two candidate models that might explain how an infectious disease 
having varying disease response could spread in a community of households. Large population 
properties of the two models are presented. These results are used to show by means of 
numerical illustrations, that it is generally possible to discriminate between the two models. 
More precisely, given final outcome data from a sufficiently large community of households it 
is, except in some degenerate cases, possible to determine which of the two explanations to 
varying disease response that best explain the data. 

Both models could of course be extended towards higher realism. For example, besides 
household structure, all individuals are assumed similar whereas it would be more realistic 
to distinguish between adults and children having different mixing rates. Another extension 
would be to allow for more than two different disease responses. It is of course also possible 
to come up with other models giving rise to mild and severe infectives. However, we believe 
that the two models studied capture the perhaps two most likely reasons: either the infection 
status of an individual is predetermined or else it depends on whom the person was infected 
by. In the first situation it could be natural to extend the model to allow this predetermined 
status to be dependent between individuals of the same household, for example due to previous 
exposure to the disease. In the present model it is assumed that the predetermined infection 
status is independent also within households. Another important extension would of course 
be to apply the method to real data with the hope to find out more about the underlying 
reason for having varying disease response. 
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